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ABSTRACT 

We investigate the influence of the ionization of helium on the low-degree acoustic oscil- 
lation frequencies in model solar-type stars. The signature in the oscillation frequencies 
characterizing the ionization-induced depression of the first adiabatic exponent 7 is 
a superposition of two decaying periodic functions of frequency v, with 'frequencies' 
that are approximately twice the acoustic depths of the centres of the He 1 and He 11 
ionization regions. That variation is probably best exhibited in the second frequency 
difference /S.2V n ,i = ^n-i,l ~ 2z/ n> ; + v n +i.i- We show how an analytic approximation 
to the variation of 7 leads to a simple representation of this oscillatory contribution 
to IS.2V which can be used to characterize the 7 variation, our intention being to use 
it as a seismic diagnostic of the helium abundance of the star. We emphasize that the 
objective is to characterize 7, not merely to find a formula for that reproduces 
the data. 
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^ 1 1 INTRODUCTION 



The issue addressed in this paper concerns the direct inter- 
pretation of an asteroseismic signature of helium ionization 
in terms of an ionization-induced property of the stratifica- 
tion of the star that produces that signature. The knowledge 
so obtained is important for designing calibrations of theo- 
retical stellar models against seismic data (together with 
other astronomical data) , for the purposes of obtaining esti- 
mates of some of the gross properties of stars, such as their 
ages and their chemical compositions. 

The first approach that might come to mind is to ad- 
just apparently appropriate parameters of the stellar models 
in such a way as to minimize some measure of deviation of 
the individual theoretical oscillation eigenfrequencies from 
the observed frequencies. Without careful regard to the ex- 
tent to which those frequencies are influenced by the stel- 
lar properties in mind, namely the helium abundance Y in 
this discussion, this naive procedure can hardly be optimal, 
although it has been used: it was first carried out for the 
Sun, by Christensen-Dalsgaard & Gough (1981), who sim- 
ply minimized the frequency misfit by least squares. The 
procedure has the merit of being simple to execute, but it 
is subject to all the uncertainties in stellar modelling, some 
of which need not necessarily be very pertinent to the ob- 
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jective of the calibration. It is more prudent to aim to work 
instead with certain combinations of frequencies that are 
insensitive to the irrelevant properties of the stellar mod- 
elling. Commonly used have been the so-called large and 
small frequency separations, which are relatively insensi- 
tive to the uncertain state of the outer layers of the stars 
(e.g. Christensen-Dalsgaard & Gough 1984; Gough 1986b). 
Christensen-Dalsgaard (1986, see also Ulrich 1986; Gough 
1987) illustrated how from a knowledge of these separa- 
tions one could in principle estimate the mass and age of 
a main-sequence star, provided that its chemical composi- 
tion were known; indeed, in the early days of helioseismology 
these separations were often preferred to the raw frequencies 
for solar calibrations (Christensen-Dalsgaard & Gough 1980; 
Shibahashi, Noels & Gabriel 1983; Ulrich & Rhodes Jr 1983). 
But the chemical composition of a star, particularly the he- 
lium abundance, is often not well known, and it is therefore 
incumbent upon us to devise more reliable ways of deter- 
mining it, especially in view of the high- precision data an- 
ticipated from the imminent space missions COROT (Con- 
vection Rotation and Planetary Transits; Baglin 2003) and 
Kepler (Borucki et al. 2003; Basri, Borucki & Koch 2005), 
and from the ground-based campaigns such as SONG (Grun- 
dahl et al. 2006) and those of the kind organized by Kjeld- 
sen, Bedding and their colleagues (e.g. Kjeldsen et al. 2005). 
Even in the Sun, in which helium was discovered by J.N. 
Lockyer and P.J.C. Janssen in prominences observed during 
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the eclipse of 1868, the spectroscopic determination of Y re- 
mains elusive to this day. This comes about because helium 
spectrum lines are formed well above the photosphere, in 
the much hotter chromosphere and corona where the devi- 
ation from local thermodynamic equilibrium can be severe, 
and where the helium abundance is not even necessarily the 
same as in the layers beneath the visible surface. 

In principle, a direct seismological signature of Y in the 
envelope of a (sufficiently cool) star can be sought through 
its effect on the sound speed, via the depression of the value 
of the first adiabatic exponent 7 induced by ionization, par- 
ticularly if the depression occurs in an essentially adiabat- 
ically stratified region of a convection zone. The variation 
is seismically abrupt, in the sense that, at least for lower- 
frequency modes, its radial extent is less than or comparable 
with the inverse radial wavenumber of the eigenfunctions. 
We call such a feature an acoustic glitch. It causes a small 
shift in the eigenfrequencies, relative to those in a corre- 
sponding putative star with no ionization, that shift being 
an oscillatory function of the frequency itself. The amplitude 
of that function increases monotonically with Y, and should 
therefore be a robust indicator of the value of Y; the fre- 
quency of the function is determined by the acoustic depth 
of the centre of ionization beneath the seismic surface of the 
star. 

Oscillatory components are produced also by other 
acoustic glitches, such as that at the base of the convection 
zone which is caused by a near-discontinuity in the gradient 
of the density scale height. They are evident in the large 
separations Ai^ nj ; = v n> i — v n -ij between the cyclic eigen- 
frequencies of seismic oscillations of order n and degree 
I, complicating the measurement of the underlying smooth 
component of Aiv n ,i, and also of other diagnostics such as 
the small separation v n .i — v n -i,i+2- Not taking them into 
account leaves an undulation in the outcome of data-fitting 
when plotted against the extremities of the range of frequen- 
cies employed (Gough 2001), which could lead to misinter- 
pretation. Therefore it is prudent to separate the smooth 
and the oscillatory components, using the two of them as 
complementary diagnostics. We address in this paper a way 
in which that separation could be carried out, by analysing 
signatures suggested by asymptotic analysis of the eigen- 
frequencies of a grid of solar models, with the intention of 
determining parameters that are good measures of the vari- 
ation of the stratification in the second helium ionization 
zone that are not unduly contaminated by other properties. 

Direct seismological signatures of the helium abun- 
dance in the solar interior have been constructed from 
intermediate-degree acoustic modes by various techniques 
(e.g. Gough 1984a; Dappenetal. 1991; Kosovichev et al. 
1992; Vorontsov, Baturin & Pamyatnykh 1992; Perez 
Hernandez & Christensen-Dalsgaard 1994b; Basu & Antia 
1995). In distant stars, however, only low-degree modes can 
be observed, and we need to find a diagnostic in only those 
modes, a need which has been stressed before by e.g. Gough 
(1990, 1998). One procedure is to attempt a direct inver- 
sion of low-degree modes for the helium abundance Y, using 
the constraint that the bulk of the convection zone is adia- 
batically stratified (Kosovichev 1993). More commonly, the 
outer phase function a(cu) in the asymptotic relation 

(n + &)tv/lo = F(w) + ... (1) 



(Gough 1984b, 1993; Vorontsov 1988; Vorontsov, Baturin 
& Pamyatnykh 1991) has been used, for it exhibits ex- 
plicitly the oscillatory components of lu induced by ion- 
ization (e.g. Baturin & Mironova 1990a,b; Christensen- 
Dalsgaard & Perez Hernandez 1992; Roxburgh & Vorontsov 
1994); here to is the angular frequency of a mode of or- 
der n and degree I, w — + \) and F is a func- 
tion of w which depends also on the structure of the star. 
The idea is to determine a by fitting the functional form 
of its asymptotic relation to the seismic data (cf. Gough 
1986a; Vorontsov 1988; Gough & Vorontsov 1995). It is ev- 
ident from that relation that a/co (and ty^ 1 F) can be de- 
termined from data u) n ,l only up to an unknown additive 
constant, so Brodsky & Vorontsov (1987) proposed using 
the derivative (3 = d(co~ 1 a)/duj, which contains the same 
information and is more elegant. The oscillatory compo- 
nent can be largely separated by appropriate filtering (cf. 
Perez Hernandez & Christensen-Dalsgaard 1994a), and the 
contribution from helium ionization used to calibrate stellar 
models to determine Y (e.g. Perez Hernandez & Christensen- 
Dalsgaard 1994b, 1998; Lopes etal. 1997; Lopes & Gough 
2001). 

An alternative procedure is to consider the second fre- 
quency difference, A2^ n ,i (Gough 1990), or even higher- 
order differences (Basu, Antia & Narasimha 1994; Basu 
1997), which are much simpler to extract from the data 
than the phase functions a and j3 yet appear to contain 
the same information. It has been found expedient to model 
their oscillatory components with parametrized functions 
derived from representations of the acoustic glitches, and 
then to determine the sought-for properties of the star by 
calibrating the parameters corresponding to a sequence of 
stellar models against the seismic data. Various approxi- 
mate formulae for the seismic signatures that are associ- 
ated with the helium ionization have been suggested and 
used, by Monteiro & Thompson (1998, 2005), Gough (2002), 
Miglio et al. (2003), Basu etal. (2004), Verner, Chaplin & 
Elsworth (2004) and Ballot, Turck-Chieze & Garcia (2004), 
not all of which are derived directly from explicit acous- 
tic glitches. Gough used an analytic function for mod- 
elling the dip in the first adiabatic exponent. In contrast, 
Monteiro & Thompson (1998) assumed a triangular form, 
which was adopted also by Monteiro & Thompson (2005) , 
Miglio et al. and Verner et al. These are the only two at- 
tempts to relate the low-degree seismic signature directly 
to the properties of the ionization zone, although the lat- 
ter appears to be too crude to capture the physics ad- 
equately (Monteiro & Thompson 2005). Basu etal. (2004), 
Ballot et al. (2004) and Piau, Ballot & Turck-Chieze (2005) 
have adopted a seismic signature for helium ionization that 
is similar to that arising from a single discontinuity in the 
sound speed (see also Basu, Antia & Narasimha 1994; Basu 
1997), which is even less appropriate; the artificial disconti- 
nuities in the sound speed and its derivatives that this and 
the triangular representations possess cause the amplitude of 
the oscillatory signal to decay with frequency too gradually, 
although that deficiency may not be immediately notice- 
able within the limited frequency range in which adequate 
asteroseismic data are or will imminently be available. How- 
ever, any inappropriate representation of the seismic signa- 
ture can hardly be trusted to measure the characteristic of 
the glitch that causes it. 
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Although the seismic quantities a, associated with 
low-degree modes and Akv n ,i have been analysed in such 
a way as to enable the magnitude of the oscillatory compo- 
nents produced by helium ionization to be obtained, in none 
of the studies, save those of Monteiro & Thompson (1998, 
2005) and Gough (2002), was the signature related directly 
to the form of the acoustic glitch that produced them. They 
were simply used as a medium for fitting theoretical models 
to seismic data, in the hope that the stellar properties of 
interest, such as Y, had been determined reliably. It has not 
been demonstrated whether the inferences are essentially un- 
contaminated by other properties of the star. Our aim in this 
paper is to design a seismic signature that truly reflects the 
properties of the seismic glitches we wish to measure. And wc 
report below that we have been reasonably successful. The 
next step in advancing our understanding would require a 
further independent study of stellar structure (and the equa- 
tion of state) to ascertain how the glitch properties relate to 
Y. The glitch amplitude increases with increasing Y . How- 
ever, we have found that the rate at which it does so depends 
also on other factors, such as the entropy of the adiabatically 
stratified ionization zones, as has been noted previously by 
Dappen, Gough & Thompson (1988), Lopes etal. (1997) and 
Perez Hernandez & Christensen-Dalsgaard (1998), and that 
therefore the glitch amplitude alone is not a simple measure 
of Y. 

In Sun-like stars, the stratification in the helium ion- 
ization zones is very close to being adiabatic (and Reynolds 
stresses are negligible), and consequently changes £7 in 7 
produce well defined changes in the variation of sound speed 
with radius. We note, in passing, that the ionization of hy- 
drogen produces much greater 7 variation, but unfortunately 
(for this study) the hydrogen ionization zone is not even ap- 
proximately adiabatically stratified throughout; the form of 
the seismic glitch must be susceptible to the uncertain for- 
malism adopted to model the convective energy and momen- 
tum fluxes, and relating it to 7 must therefore be unreliable. 
For that reason we concentrate on the ionization of helium, 
as have others before us. We represent the variation of 7 
with a pair of Gaussian functions. This correctly results in a 
decay of the amplitude of the seismic signature with oscilla- 
tion frequency that is faster than that which the triangular 
and the single-discontinuity approximations imply, and also 
takes some account of the two ionization states of helium. 

The plan of the paper is as follows: in the following Sec- 
tion we introduce and derive a simple seismic diagnostic of 
a Gaussian acoustic glitch produced by only the (dominant) 
second ionization of helium, together with a representation 
of the contribution from the base of the convection zone, 
the latter being derived in Appendix B; the variation of the 
eigenfunction phase tp with acoustic depth r is first repre- 
sented by the simple formula tp ~ uit + e, where e is taken 
to be constant. This diagnostic is tested against artificial 
frequency data computed from two sequences of theoretical 
solar models, which we describe in Section 3. Although the 
diagnostic can be fitted to the theoretical eigenfrequencies 
tolerably well, the inferred value of the acoustic depth m of 
the glitch is too low, and the inferred amplitude — 8^/~i\ T=Tll 
is too high. This is as one might expect, partly because in 
this initial phase of the investigation we have adopted only 
a single Gaussian function to represent the two ionization 
stages of helium, and partly because we have overestimated 



-!/>(th) by not taking account of the acoustic cutoff frequency 
w a . We rectify this, in two steps, in Section 4, first by incor- 
porating w a into V>, which essentially halves the discrepancy 
between the true Tn and its seismologically inferred value, 
and also reduces the amplitude discrepancy considerably (al- 
though it does not reduce the measure \ 2 of the data misfit), 
and then by adding in a simple way a contribution to the 
diagnostic from the first stage of ionization of helium, which 
reduces the An discrepancy by a factor three (and also re- 
duces x )• We must point out, however, that the value we 
must adopt for the ratio Ai/An of the characteristic widths 
of the He I and He II glitches in order to obtain a good fit to 
the artificial data is rather larger than that suggested by the 
models. We do not know why, although a change of this kind 
is perhaps suggested by the functional form of the contribu- 
tion to the depression of 7 due to the ionization of hydrogen, 
which is depicted in Fig. 7. We discuss our analyses further 
in Section 5, and we draw our conclusions in Section 6. 



2 THE SEISMIC DIAGNOSTIC 

A convenient and easily evaluated measure of the oscillatory 
component is the second difference with respect to order n 
of the cyclic frequencies v n ,l of like degree I: 

A 2 ^„,i = Vn-1,1 - 2Vn,l + Vn + 1,1 , (2) 

(Gough 1990). This measure is contaminated less than the 
first difference Aiz/ nj ; = v n ,i — v n -i,l, otherwise known as the 
large frequency separation, by the smoothly varying compo- 
nents of A\v (here and henceforth we simplify the notation, 
where it is not ambiguous to do so, by omitting the sub- 
scripts n, I). And it is less susceptible to data errors than the 
fourth and other higher-order differences (see Appendix C); 
moreover, adopting higher-order differences risks requiring 
more consecutive modes than might be available to provide 
a practical diagnostic. 

Any localized region of rapid variation of the propaga- 
tion speed of an acoustic wave, which here we call an acous- 
tic glitch, induces an oscillatory component in A 2 ^ with a 
'cyclic frequency' approximately equal to twice the acoustic 
depth 

r R 

t= c _1 dr (3) 

*' r glitch 

of the glitch, where r is a radial coordinate, c is the adia- 
batic sound speed and R is the radius of the seismic surface 
of the star (essentially the surface on which c 2 , if extrapo- 
lated linearly outwards from the outer layers of the adiabat- 
ically stratified region of the convection zone, would vanish 
- cf. Balmforth & Gough 1990). The amplitude of the oscil- 
latory component depends on the amplitude of the glitch, 
and decays with v once the inverse radial wavenumber of the 
mode becomes comparable with or greater than the radial 
extent of the glitch. By calibrating a theoretical representa- 
tion of the effect of glitches against the observations one can, 
in principle, learn about the characteristics of those glitches. 
Because only low-degree modes are used, the /-dependence 
can safely be ignored, even down to the base of the convec- 
tion zone (except, possibly, in very- low-mass, almost fully 
convective stars). We remark further on that in Section 5. 
If a frame of reference exists in which the basic structure of 
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a star, which in general is presumed to be rotating (but only 
slowly) , is independent of time, then one can seek linearized 
adiabatic oscillations in that frame whose time dependence 
is sinusoidal with angular frequency oj. The adiabatic oscil- 
lation frequency u> is related to the displacement eigenfunc- 
tion £ through a variational principle (Chandrasekhar 1963; 
Lynden-Bell & Ostriker 1967), the derivation of the acous- 
tic signature from which we now address. We consider the 
structure of the star to be described by smoothly varying 
functions upon which are superposed glitches, each varying 
on a small spatial scale. The contribution Slo from those 
glitches to the eigenfrequency u of a mode of oscillation 
(with vanishing pressure gradient at the surface), can then 
be written (see Appendix A) 

. SIC — ujST , . 

(4) 

where X is the mode inertia and SIC depends on the dis- 
placement eigenfunction £ and on the glitches in the equi- 
librium quantities, such as the (localized) perturbation S'y 
to the first adiabatic exponent 7 = (d lnp/dln p) a , s be- 
ing specific entropy, caused by ionization. In equation (4), 
contributions from the surface (which at most contribute a 
relatively smoothly varying component to A2V) have been 
ignored. 

We separate Slj into a smoothly varying component 
8sm<jJ and an oscillatory component 5 obc uj: 



Suj = S sm LU + S 



(5) 



(and we adopt a similar notation for other variables) and 
we approximate, with a suitable eigenfunction normaliza- 
tion, the oscillatory component 8 ltOSC u} associated with he- 
lium ionization as follows: 



c u> — osc(S-yLu) , 



where 



5-yU) 



SjlC 



with 

SylC ~ SjICi 



(Sy) p(div£) 2 r 2 dr , 



(cf. equations (A5),(A6)) and 



I • I r dr , 



(6) 



(7) 



(8) 



the integrations being over the entire extent of the star, 
where p and p are pressure and density of the equilibrium 
state, and S'y is a suitably defined rapidly varying perturba- 
tion to 7 induced by ionization; osc denotes oscillatory part. 
We have confirmed numerically that the terms neglected in 
SIC are indeed substantially smaller than S^IC (although they 
may not be wholely negligible, as is evinced by the differ- 
ences in 57/7 1 t u in Figs 9 and 10). In view of the varia- 
tional principle, we have approximated the eigenfrequencies 
by those of a corresponding nonrotating smoothly varying 
(i.e. having no glitches) model star; consequently in equa- 
tions (7) and (8) we have been able to adopt the usual sepa- 
ration of the components of the eigenfunctions £ in terms of 
functions of r and spherical harmonics, carry out the angular 
integrations, and regard £ to be a function of r alone. In the 
asymptotic limit of high-order modes the (suitably normal- 
ized) radial component £ of the displacement eigenfunction 



£ and the associated Lagrangian pressure perturbation p are 
given by (e.g. Gough 1993) 



rp 



1/2 



u;V /2 



where the phase function ip is given by 



K(r') dr' + -~wr + e, 



(9) 



(10) 



in which here we have approximated the radial component 
of the wavenumber of high-order modes as K ~ iv/c, and 
we have assumed that r is not close to rt, the upper turning 
point at which the mode is reflected; e is a phase constant 
(strictly speaking, a phase variable, for it varies slowly with 
ui, but in the simple approach that we adopt in just this 
section that variation is ignored) which takes some account 
of the fact that K deviates substantially from uj/c near the 
upper turning point. It follows that 



(div|) 2 = ( ^) ~ - sin 2 ip , 

and therefore equation (7) can be written in terms of 
<?7 



^7,osc/Cl 



1 f „ 

— -us I — cos 2(ujt + e) dr , 

2 J 7 



(11) 



(12) 



in which the limits of integration are such as to include the 
region in which £7 is nonzero; moreover 



f T 1 
T ~ uj / cos (ujt + e) dr ~ —Tu , 

Jo 2 



(13) 



where T is the acoustic radius of the star (i.e. the sound 
travel time from the acoustic surface to the centre of the 
star). Note that the normalization of this expression for the 
inertia differs from the usual. The neglect of the acoustic 
cutoff frequency in the expression for K will be dropped 
in Section 4.1. However, we shall continue to neglect the 
degree-dependence, which is a good approximation for the 
low-degree modes detectable by asteroseismology; we justify 
that neglect in Section 5. 

In this introductory discussion we represent the acous- 
tic glitch induced by the (second) ionization of helium by 
a Gaussian function about the acoustic depth r = rn of 
the centre of the He II ionization region (see Fig. 3) , as did 
Gough (2002): 



^ ~G(r;r„,r„,A„) 

7 



j_rn 

2T Ai 



-(t-t„) 2 /2A„ 2 



(14) 



in which Tu and An are constants and An is much smaller 
than both m and T — m. We improve on this simple ex- 
pression in Section 4.2. The oscillatory component, 5 7iOSC u;, 
imparted by the glitch can then be estimated asymptotically 
from equations (6), (12) and (13) to be 



5 JtOBC u) ~ Anu e 2An w cos 2(Tnw + en) • 



(15) 



in which we have replaced the phase e by en, and we have 
introduced a normalized glitch amplitude 



An = lr u T 



(16) 



Notice that the amplitude of the oscillatory component 

which for large ui 2 



5 7i osc^ decays exponentially with ui 2 
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1000 1500 2000 2500 3000 3500 4000 
v (/xHz) 

Figure 1. The symbols are second differences A2^, defined by 
equation (2), of low-degree solar frequencies obtained from GOLF 
data and kindly supplied by J. Ballot. Most of those data have 
been presented by Ballot et al. (2004). The most extreme outly- 
ing datum has been removed. The vertical bars represent stan- 
dard data errors, evaluated under the assumption that the er- 
rors in the raw frequencies are independent; their effective over- 
all value is (cr) = 4.3 nHz. The curve is the seismic diagnos- 
tic Do(y;ak), which has been fitted to the data in a manner 
intended to provide an optimal estimate of the eleven param- 
eters a k . The values of some of those fitting parameters are: 
-hh\r u = Ah/v^i/qAii ~ 0.065, m ^ 707s, A n =± 52 s. 
The measure E of the overall misfit is 50 nHz. The direct mea- 
sure x is 9.9; the minimum-x 2 fit of the function Do to the data 



yields x n 



7.0. 



is faster than the algebraic decay predicted by the repre- 
sentations of glitches containing discontinuities (cf. equa- 
tion (17)); the amplitudes due to those perturbations decay 
only as u)~ q , where q is the order of the lowest derivative of 
a dynamically relevant variable of the equilibrium state that 
contains a discontinuity. 

The acoustic glitch at the base of the convection zone 
(r = r c ) also contributes an oscillatory component, S CtOBC uj, 
to the frequency ui of the modes, and hence to A2V. This 
component must be included when fitting functional forms 
for A2V to real data, even though, in this investigation, our 
interest is restricted to helium ionization. The glitch is es- 
sentially a discontinuity in the acoustic cutoff frequency, and 
is reflected as a near discontinuity in d£/dr, but not in £ 
(a local mixing-length model - as we have adopted in the 
construction of our theoretical test models - leads to a true 
discontinuity) . From Appendix B we derive the following ap- 
proximate expression (see equation (B19)); 



with 



, 3 — 2 / -1 1 -1/4 2 2\ -1/2 

i AcLOqUJ (1 + 1/4t lu j 
x cos [2(t c o> + e c ) + tan -1 (2too>)] , 



d 2 lnp 



dr 2 



(17) 



(18) 



where r c is the acoustic depth of the base of the convection 
zone and c c = c(r c ); To <§C T is a measure of the character- 
istic acoustic distance in the radiative interior over which 
the glitch causes the acoustic cutoff frequency to deviate 
substantially from a smooth function, and was determined 



to be 80s by fitting an exponential function to a theoreti- 
cal reference solar model (see Appendix B). The (constant) 
quantity uio is the asymptotic mean large frequency separa- 
tion: 



my 



dr 



7T 

T 



(19) 



which was determined from fitting by least-squares the 
asymptotic expression for the oscillation frequencies (e.g. 
Tassoul 1980; Gough 1986b) to the frequency data, namely 



~ (n + § I + s)vo - 



Bl{l + \)-C 2 



v n ,l 



(20) 



in which vq — loqI2-k and e, B and C are also constants. 
The amplitude A c is a measure of the discontinuity in the 
second derivative of p at the base of the convection zone. The 
amplitude factor (1 + is only a weak function 

of w, especially above cyclic frequencies of about 3 mHz in 
the solar case. 

The amplitudes An and A c , and the phase increments 
2en and 2e c + tan -1 (2tooi), vary only slowly with u>. This 
property facilitates the evaluation of the second differences 
of the oscillatory frequency perturbations given by equations 
(15) and (17), which is described in AppendixC, in which 
en and e c are taken to be constants. Those constants can- 
not be taken to be the same because between m and r c the 
relation K — lu/c is not exact. Indeed, this property also 
provides part of the reason why these constants cannot be 
written in terms of the phase e appearing in the asymptotic 
eigenfrequency equation (20). Smooth contributions to the 
second differences of the eigenfrequencies (arising, in part, 
from refraction in the stellar core, from hydrogen ionization 
and from the superadiabaticity of the upper boundary layer 
of the convection zone which lies in an evanescent region of 
the acoustic mode) must also be accounted for. Guided by 
the functional form of the asymptotic relation (20), we ap- 
proximate them by a third-degree polynomial in lj^ 1 . They 
account for both the second frequency differences of a pu- 
tative reference solar model that has no glitches associated 
with helium ionization or the base of the convection zone, 
and the smooth contributions from those glitches. The com- 
plete fitting function is then given by 



A 2 v ~ A 2 (<5 7 , osc ^ + 5 C 



fc=0 



(21) 



i FyiAyiv c~ S7r AuV cos[2(27rni^ + en) - Sn] 

+ F C A C ^- 2 (1 + l/167r 2 roV)~ 1/2 

x cos[2(27tt c !/ + e c ) + tan _1 (47rro!/) — S c ] 

3 

+ ^2 a k v~ k = D {v; a k ) , 



(22) 



where we have now converted to cyclic frequency v, which 
we regard as a continuous variable; formulae for the ampli- 
tude factors Fu and F c and for the phase increments 5n and 
<5 C arising from taking the second differences are given by 
equations (C6) and (C7) with a and b given immediately 
following equation (C9) in the case of Fn and Su and by 
equations (C12) in the case of F c and <5 C . The quantities 
a fe are the eleven parameters An, An, m, en, A c , r c , e c and 
a k (k = 0, 3). 
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Figure 2. Denotation of the nine calibrated solar models which 
are used for testing the formulation of the second frequency differ- 
ences. The 'central model' is model 0; the set of the four models 
(1-4) have a constant value of Z but varying age; the second set 
of the models (5-8) have constant age but varying Z. 



In Fig. 1, second differences A2V of solar oscillation 
frequencies obtained from the Global Oscillations at Low 
Frequency (GOLF) investigation on the Solar and Helio- 
spheric Observatory (SOHO) spacecraft are plotted against 
v. For comparison, the seismic diagnostic Do is plotted 
as a continuous curve against v, which has been fitted to 
the data with the intention of providing an optimal es- 
timate of the eleven parameters otk by minimizing E 2 — 
{Do — d) t C 1 {Do — d) I 1 , where d is the vector whose 

components are the N data iX%Vn,i, Do is the iV-dimensional 
vector with components Do(y n ,i'> a k)> C is the covariance 
matrix of the errors in the data d, and the superscript t de- 
notes transpose. The covariance matrix C, whose eigenvalues 
are denoted by A*,, was evaluated under the assumption that 
the errors in the measurements of the frequencies v n j are in- 
dependent random variables with variance cr^ we take the 
square root of the harmonic mean of At, as an overall measure 
(a) of the magnitude of the errors in the second differences. 
A direct measure of the fit of the curve to the data is given 
by X 2 = N- 1 E„,i{[- D o(^„, i ;Qfc) - A 2 v„j]/a n ,i} 2 - 

A similar comparison may be made with comparable 
second differences of frequencies obtained from the Michel- 
son Doppler Imager (MDI), also on SOHO, which, over 
the available frequency range, 1473-3711 /^Hz, of those data 
(from which, for equity, we have also removed the most ex- 
treme outlier) have overall error (a) = 5.3 nHz, and deviate 
from the fitted diagnostic (22) by E — 32 nHz, compared 
with (a) = 6.4 nHz and E = 34 nHz for the GOLF data 
over the same frequency range. The sharp upturn of the 
purported 'smooth' contribution (dashed curve) in Fig. 1 as 
v decreases below 1300 /iHz is probably an artefact of the 
failure of the asymptotic approximation upon which equa- 
tion (22) is founded; this deficiency is removed partially by 
the improved diagnostic D2 which we introduce in Section 4. 

One should expect equation (22) to represent the sec- 
ond differences of the actual frequencies more faithfully 
than the expressions that have been used before (e.g., Mon- 
teiro & Thompson 1998, 2005; Gough 2002; Basuetal. 2004). 
This is so partly because here we use a representation of the 
Hell glitch in 7 that is more realistic than the non-analytic 
functions adopted by Monteiro et al. and Basuetal., and we 
also use a more realistic representation of the stratification 
immediately beneath the convection zone. 



1.64 

1.62 

1.60 

1.58 

1.56 

1.54 

1.52 
1.50 



7n_ 



700 



800 



900 



1000 



T ( S ) 



Figure 3. Adiabatic exponent 7 (solid and dotted curve) for 
the central stellar model as a function of the acoustic depth 
t, which is measured from the acoustic surface of the model 
(estimated as the place where the linear extrapolation of c 2 
from the region of the upper turning points of the modes van- 
ishes, some 225 s above the photosphere). The dashed curve is 
7110 - (7110 - 7ll) ex P[-( T - T ll) 2 /2A 2 I ], where 700 = 1.651 is 
the value of 7(p(tji), i9(tti), 0), evaluated at Y = (dot-dashed 
curve in Fig. 7) and is temperature; the remaining parameters 
(7n, Tu, ^II) have been adjusted to fit 7 in the region in which 
the curve representing 7 is solid, namely m — Tf < T < Trr -f Tf 
with Tf ~ 30 s. The parameter 7jj is the minimum value of 7 in 
the He II ionization zone, and ttt is the acoustic depth at that 
minimum. The parameters associated with the other test models 
were determined likewise, with values of Tf as close as possible to 
that of the central model subject to their fitting exactly onto the 
computational mesh; their values differ from that of the central 
model by less than 3 s. 



3 TESTING THE FORMULATION ON 
THEORETICAL MODELS 

3.1 Test models 

Two sets of five calibrated evolutionary models for the 
Sun were used. The models had been computed by 
Gough & Novotny (1990) using the procedure described by 
Christensen-Dalsgaard (1981). The opacity values were ob- 
tained from interpolating in the Los Alamos Opacity Li- 
brary for the Grevesse (1984) mixture, supplemented at low 
temperature by interpolating in the tables of Cox & Steward 
(1970a,b). Gravitational settling was not included. One set 
of models has a constant value for the heavy-element abun- 
dance Z but varying age in (logarithmically) uniform steps 
of magnitude approximately 0.051 in \nt (t being age); the 
other has constant age but varying Z, in uniform steps of 
magnitude approximately 0.016 in In Z. The central mod- 
els (model 0, see Fig. 2) of the two sets coincide, having 
t = 4.60 Gy and Z = 0.02014. All models were calibrated by 
adjusting the initial helium abundance Yo and the mixing- 
length parameter to satisfy the solar values of the luminos- 
ity and photospheric radius: Lq = 3.845 x 10 33 ergs _1 and 
R Q = 6.9626 x 10 10 cm, resulting in a central model with 
Yo = 0.2762. The models were examined carefully and ad- 
justed to satisfy hydrostatic equilibrium to high precision. 
Eigenfrequencies of adiabatic oscillation modes with I < 2 
were computed for all models in the Cowling approximation 
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(gravitational forces are long-range and do not materially 
affect glitch diagnostics, so the Cowling approximation is 
quite adequate for our purposes) using the Aarhus Adiabatic 
Pulsation Package 1 . From these were constructed second fre- 
quency differences according to equation (2). 

3.2 Test results 

For all nine test models we determined the properties of the 
He II depression in 7 by fitting a Gaussian (plus a constant) 
to 7. Details of this fitting process are illustrated in Fig. 3 
for the central model 0. The three parameters, — 5^//^/]^ , Tn 
and An so determined are plotted in the top three panels of 
Fig. 4 for the nine test models as a function of the helium 
abundance Y, where -Sj/j\ Tll = 2(7110 - 7ii)/(7iio + 711). 
The value of 7110 is the value of 7 evaluated with the pressure 
and temperature of the model at r = Tn, but with a helium 
abundance Y = (see dotted curve in Fig. 7). For the central 
model 7110 = 1.651. The values of the other models differ 
from that of the central model by less than 2xlCP 4 . 

The first two lower panels in Fig. 4 display the value of 
A c (equation (18)) and the acoustic depth r c of the base of 
the convection zone. These values were obtained by fitting 
analytical functions (as discussed by Christensen-Dalsgaard, 
Gough & Thompson 1991) to the first logarithmic density 
derivative, dlnp/dlnr, of the models either side of the dis- 
continuity at the base of the convection zone, as explained 
in Appendix D. The third (rightmost) panel is a standard 
measure (mean squared residuals) \ 2 f° r the goodness of 
the Gaussian fit to the Hell depression in 7. 

The asymptotic formula Do was then fitted to the sec- 
ond differences of the computed eigenfrequencies of the mod- 
els as described in Section 2; some of the resulting param- 
eters are plotted in Fig. 5. There is a superficial resem- 
blance to the corresponding direct estimates from the mod- 
els depicted in Fig. 4 (except for the slope of Au(Y) for 
the Z- varying sequence, which we suspect is influenced sub- 
stantially by heavy-element ionization, and, of course, x 2 > 
which has quite different meanings in the two figures) . Per- 
haps the most obvious other difference is that the slopes of 
the S~//j\ TlI curves from the two sequences of direct model 
measurements are almost the same, whereas those of the 
frequency- fitting are not. The magnitudes of the values of 
some of the parameters are also discrepant. As will become 
evident in Section 4, that is to be expected in the case of 
the parameters characterizing the ionization, partly because 
the seismological fit is influenced by the effect of He I ion- 
ization, which has been treated erroneously as being part 
of He II. The acoustic depths of the bases of the convection 
zones agree reasonably well. But the amplitudes do not; we 
do not know why. 



4 FURTHER IMPROVEMENTS 

4.1 Improved treatment of the asymptotic phase 
function ij> 

One of the obvious deficiencies of our introductory analysis 
described in the previous section is our treatment of the 

1 http://astro.phys.au.dk/~jcd/adipack.11 



phase function ip, which we approximated by tp ~ lot + e, 
neglecting the acoustic cutoff frequency tj a (and also the 
relatively small horizontal component of the wavenumber) 
in evaluating the vertical component K of the wavenumber, 
whose full planar form is given by equation (B2). In this 
section we account for the effect of the frequency dependence 
of the upper turning point by keeping w a : 

representing the stratification of the outer stellar layers lo- 
cally by a polytrope of index m, for which the acoustic cutoff 
frequency associated with the Lagrangian pressure pertur- 
bation (and corrected for the singularity at the top of the 
polytropic envelope - see Gough 1993) is 

u; a = (rn + l)/r . (24) 

Thus, recognizing that the acoustic glitch associated with 
the helium ionization is confined to a small range of acoustic 
depth, we expand the phase function V>(t) about the (phase- 
shifted) centre of the glitch, r = Tn, having accommodated 
the phase constant en by replacing loth by loth = loth + en, 
yielding 

V>(t) ~ V(tii) + lok u (t - rn) + 0[(r - fn) 2 ] , (25) 
with Kn = k(tii), where 

We adopt m = 3.5 for all models; it was obtained from 
the slope of I/lo^t) in the outer layers of the adiabatically 
stratified region of the convection zone in the central model, 
and relating that to m using equation (24). We note, in pass- 
ing, that this value is close to 3.0, the effective value of m 
in the region of the upper turning points of solar modes of 
intermediate degree determined by calibrating the asymp- 
totic eigenfrequency formula (1) against observation (e.g. 
Christensen-Dalsgaard et al. 1985); actually, the values of x 2 
of the data-fits that follow are reduced somewhat if the value 
3.5 is replaced by 3.0, but then the overestimates of cor- 
responding values of the amplitudes — <57/7| T „ are greater. 
Taking expansion (25) to higher order makes a correction of 
magnitude typically less than one per cent, so we keep only 
these first two terms. The leading term in the expansion (25) 
is given by 

r(m+l)/w n 

4>(f u ) ~ lo kAt + - 

J Til 

, ,n _i / m + 1 \ -K ,„_. 
= kiitiiw - (m + 1) cos — I + — ■ (27) 

V TULO J 4 

The analysis proceeds as in Section 2 and Appendix A. The 
outcome is to change the inertia X to 

X ~ 1>{T)-j ~ iTc-i(m + l)7r 

- ^[" + 3( m + 1 H" 1 ■ ( 28 > 

The contribution to v from the glitch then becomes 

V ^ Anitf [v + ±(m + l) I / ][e" 8 ' r2 ' t " A "'' 2 cos(2V'ii)-l] , 

(29) 
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Figure 4. The top panels show — <57/7|r n := 2(7no — 7ii)/(7iio + 7n)> T II an d An, which were determined from fitting a Gaussian 
(sec Fig. 3) to the He II depression in 7, for all nine test models, plotted against helium abundance Y. Models with varying age and 
constant heavy element abundance Z are connected with solid lines; models with varying Z and constant age are connected with dashed 
lines. The central model and the models 1, 4, 5 and 8 are indicated in the upper left panel. Also in the upper left panel is a dotted 
straight line from the origin through the central model; it is evident that the effective amplitude of the helium-induced acoustic glitch is 
consistent with it being a linear function of Y in the Z-varying sequence, but not strictly so in the sequence with varying age. The lower 
right panel is the standard measure \ 2 (normalized integral of squared residuals: \ 2 = (2Ajj) 1 / (<57/7 — G) 2 dr, in which the integral 
is evaluated from Tn — An to tjj + An) of the goodness of the fit between 7 and the calibrated Gaussian. The first two lower panels 
show the properties of the density discontinuity (amplitude A c and acoustic depth r c ) at the base of the convection zone, which were 
determined from the equilibrium models. 
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Figure 5. Model properties —5"i/i\tu = An/v^rfo An, Tn, An, A c and t c for all nine test models determined from the diagnostic 
Dq defined by equation (22). The fitting parameters An, Tn, An, en, A c , t c and e c have been adjusted to fit the second differences 
A2U, defined by equation (2), of the low-degree (£=0,1,2), adiabatically computed, model eigenfrcqucncics. Note that the inferred values 
of — <57/7|tti actually decrease, gradually, with increasing Y in the Z-varying sequence. The frequency range used in the least-squares 
fitting is that available to the BiSON data, with which we compare our final diagnostic formula in Fig. 12: it is 1322-4058 /iHz. The 
lower right panel is the standard measure \ 2 (mean squared residuals, in units of (/iHz) 2 ) of the goodness of the fit between A2U of the 
modelled frequencies and the calibrated diagnostic. The line styles are as in Fig. 4. The dotted line in the upper left panel is again a 
straight line from the origin; for comparison the thin solid line is drawn with the same slope as the age-varying sequence in Fig. 4. 
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Figure 6. Model properties — S^j/^m, Tn, An, A c and r c determined from the diagnostic Di defined by equations (21), (29) & (30) 
which includes an improved treatment of the phase function ip given by equation (27). Results are shown for the nine test models. The 
fitting parameters An, Tn, An, en, A c , r c and e c have been adjusted to fit the second differences A2U, defined by equation (2), of 
the low-degree (1=0,1,2), adiabatically computed, model frequencies. The smallest frequency value used in the least-squares fitting is 
v ~ 1322 /liHz, the largest is 4058 /iHz. The lower right panel is the standard measure \ 2 (mean squared residuals, in units of (/1Hz) 2 ) 
of the goodness of the fit between A2f of the modelled frequencies and the calibrated diagnostic. The line styles are as in Fig. 4. The 
dotted line in the upper left panel is a straight line from the origin; the thin solid line is drawn with the same slope as the age-varying 
sequence in Fig. 4. 



in which An is still given by equation (16), and ipn = ip(fii), 
with V>(fn) given by equation (27). The effect on the oscil- 
latory component <5 7jOSC ^ is both to modify the amplitude 
of the cosine in equation (15) and effectively to impart an 
additional frequency-dependent contribution to the phase. 

Beneath the convection zone the cd a correction is small, 
and can be ignored; for the Sun, for example, 1 — n(r) < 
1 x 1CP 2 for r > t c . Moreover, the measure!* of the inertia, 
introduced in Appendix B, is unchanged. Consequently the 
contribution from the acoustic glitch at the base of the con- 
vection zone is influenced significantly by cd a only through 
the phase ip c = tp(f c ), where u>f c = ujt + e c , leading to 

8 c v ~ A c vlv~ 2 (l + 1/16tt 2 tqv 2 ) 1/2 

x |cos[2t/> c +tan~ 1 (47TTo^)]-(167r 2 To*/ 2 + i) 1/2 } , (30) 

where A c is given by equation (18) and tp c by equation (B20). 

The improved diagnostic function D\{y;a.k) for repre- 
senting the second frequency differences associated with the 
acoustic glitches is then given by equation (22) with the 
first term modified according to the oscillatory component 
of equation (29) and the second term by equation (30), and 
still with Fu and Su defined by equations (C6) and (C7) but 
with the extended expressions for a, a and b given by equa- 
tions (C10) and (Cll). The form of the seismic diagnostic 
from the base of the convection zone is unaffected, save for 
the replacement of 2-kt c v + e c by tp c . 

Results for the seismic diagnostic D± are shown in Fig. 6 
for all nine test models (in which we have, as before, consid- 
ered the smooth terms - i.e. the last terms in both equations 
(29) & (30) - to have been incorporated into the polynomial 



in Compared with the results of the previous section, 

plotted in Fig. 5, the amplitudes — £7/7!^ are smaller and 
the values of Tn, An and r c are larger, all being closer to 
the model values depicted in Fig. 4. 

4.2 Including He I ionization 

A second deficiency in our initial procedure is that we took 
no explicit account of the first ionization of helium. The 
reason is that He I ionization is merged with the ionization 
of hydrogen, which dominates in the depression of 7 even 
at the centre of He I ionization. However, the effect of 
the first ionization of helium is, as is evident in Fig. 7, 
not negligible; it broadens the hydrogen- induced dip in 
7. Without precise knowledge of the hydrogen-induced 
profile (perhaps even with it) seismological calibration of 
the broadening is difficult. Our approach is not even to try. 
Instead we simply attempt to relate directly the reduction 
of 7 by He I ionization to that due to Hell ionization, which 
we do aim to measure. The much greater effect of hydrogen 
ionization is, on the whole, acoustically shallower than the 
effect of He I, and therefore produces a smoother seismic 
signal, and we trust that it is adequately accommodated by 
the polynomial ikV~ k in the expression (21) for A2^ . 

The degree of helium ionization, and its effect on 7, is 
determined by the law of mass action for ionizing species, 
which is approximated by Saha-like equations (we use the 
approximate equation of state of Eggleton, Faulkner & Flan- 
nery 1973). It depends not only on the temperature and 
the number density of helium, but also on the free-electron 
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Figure 7. The solid curve is the adiabatic exponent 7(p,#, Y), 
where is temperature, through the ionization zone of the central 
model 0; the dot-dashed curve is the adiabatic exponent ~f(p, 0). 
The dashed curve is the sum of the first two terms of the Taylor 
expansion of 7 about Y = given by equation (31). Note that 
the amplitude of variation of the dashed curve in the He II ion- 
ization zone is greater than that of the solid curve; the difference 
contributes to the discrepancy between the values of S-f/^Tu in 
Figs 9 and 10. 



density, which itself depends on the abundances of all ion- 
ized species, especially hydrogen. In general these factors are 
nonlinear ly dependent on one another, but we base our pro- 
cedure here on the idea that because the number density of 
helium is quite low, although not extremely low, the influ- 
ence on 7 is approximately a linear function of the helium 
abundance Y. As we demonstrate below, the linear approx- 
imation works quite well, although it is evident from the 
top left panels of Figs 4, 5, 6, 9 and 10 that the linear ap- 
proximation is not strictly correct. It is important to realize, 
however, that it is used solely for estimating the He I con- 
tribution to 7 relative to the larger He II contribution, and 
therefore plays only a minor role in calibrating Y : our expec- 
tation is that it can account adequately for the deviation of 
the form of the seismic helium signature from those implied 
by the simple equations (15) and (29), thereby producing a 
more faithful match with the data, and consequently per- 
mitting a more robust calibration of the models against the 
star. 

We first consider 7 to be Taylor-expanded about Y = in 
a reference stellar model, yielding 



7(T,y)~ 7 (T,0) +d~f/dY\ 



Y . 



(31) 



where the coefficients 7 and dj/dY on the right-hand side of 
equation (31) are evaluated at the temperature and pressure 
of the reference model but with Y = 0. We hold tempera- 
ture, rather than density, fixed under differentiation, be- 
cause ionization depends much more sensitively on We do 
not take heavy elements explicitly into account here because 
d^/dZ\z=o, where Z is the total heavy element abundance, 
is rather smaller in magnitude than d^f /dY\y=o, and in any 
case Z itself is much smaller than Y . 

The variation £7 in the adiabatic exponent 7 induced 
by helium ionization then becomes: 



Figure 8. The solid and dotted curve is d^/dY\y = o f° r tne cen- 
tral model 0, the derivative being taken at constant p and 1? (and 
Z). The dashed curve is the sum of the two Gaussian functions 
used in the construction of the seismic diagnostic (37) whose pa- 
rameters have been adjusted to fit 8^/8Y\y = o in the region in 
which the curve is solid. 



Sy ~ dj/dY\ Y . 



(32) 



The first-order expansion (31) is illustrated in Fig. 7 for the 
central model 0; the solid curve denotes the adiabatic expo- 
nent 7(p, 1?, y) of the model, the dotted curve is y(p, i?,0), 
and the dashed curve is the right-hand side of equation (31). 
A plot of d~//dY\Y=o is presented in Fig. 8, also for the cen- 
tral model (solid and dotted curve) , which clearly exhibits 
the two glitches induced by He I and Hell ionization, and 
which we formally represent by the sum of two Gaussian 
functions about the acoustic depths r = n (He 1) and t = Tn 
(Hell): 



c?7 I 
dY\ 



2-kY 



Ti 



+ A„ e 



(r-r„) 2 /2A^ : 



(33) 



in which 7 and Y are the actual values in the original model. 
The dashed curve in Fig. 8 is a least-squares fit of the 
two Gaussian functions to d^f / dY\y=o of the central stel- 
lar model 0. 

The fitting parameters for the nine reference models are il- 
lustrated in Fig. 9. Note that the ratios (3 = FiAh/ThAi, 
77 = ti/th and fx — Ai/Au hardly vary. Indeed, they vary at 
most by about 35% amongst adiabatically stratified model 
envelopes whose masses and radii vary by factors of five. 
Therefore, for the purpose of including He I ionization, we 
regard them as being constant. We adopt the values [3 = 0.45 
and r] — 0.7. But as is the case of several of the other pa- 
rameters, particularly An, we have found that the value of 
fx that gives the best fit to the frequency data is not the 
same as the one that we have attempted to measure directly 
from the models. This is no doubt a result of contamination 
by hydrogen ionization, and perhaps ionization of heavy el- 
ements too. So instead we have adopted fx = 0.90 for all 
models, which gives an optimal fit. 

The asymptotic evaluation of the contribution to Scu 
from He I ionization proceeds analogously to the treatment 
of He II, except that now one must recognize that for low- 
frequency modes the upper turning point r = r t can be 
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Figure 9. Properties of the Hell depression in 7 for all nine test models as a function of Y. The Hell properties (top panels), including 
Tu, were obtained from fitting the two Gaussians (equation (33)) to d-f/dY\Y=o of the nine test models; they are somewhat greater 
than the acoustic depths of the minima in 7 plotted in Fig. 4, apparently because the locations of those minima are influenced by the 
T-dependent contribution from the hydrogen ionization. Partly as a consequence of fitting the Gaussians over an extended range of r, 
the values of — 57/7^ , which were computed according to equation (32), are somewhat greater than those in Fig. 4, possibly the result 
of a contribution from the ionization of the heavy elements; however, the difference is not great (cf. Fig. 7). The values (not plotted here) 
of A c and r c are, of course, the same as those in Fig. 4. The line styles are as in Fig. 4. Also in the upper left panel is a dotted straight 
line from the origin through the central model; the effective amplitude of the helium-induced acoustic glitch is consistent with it being 
a linear function of Y in the Z-varying sequence, but not strictly so in the sequence with varying age. The lower panels show the values 
of fi, i] and \i. 
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Figure 10. Model properties — S"f/j\ Tn , defined as in Fig. 5, m, An, A c and r c determined from the seismic diagnostic D2 defined 
by equations (21), (30) & (37), which includes a description for both the He I and the Hell ionization zones. Results are displayed 
for the nine test models. The fitting parameters An, in, An, til, A c , t c and e c have been adjusted to fit the second differences A2^, 
defined by equation (2), of the low-degree (£=0,1,2), adiabatically computed, model frequencies. The smallest frequency value used in 
the least-squares fitting is u ~ 1322 /iHz, the largest is 4058 ^iHz. The lower right panel is the standard measure \ 2 of the goodness of 
the fit between A2V of the modelled frequencies and the calibrated seismic diagnostic. The line styles are as in Fig. 4. The dotted line 
in the upper left panel is a straight line from the origin; the thin solid line is drawn with the same slope as the age-varying sequence in 
Fig. 9. 
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Figure 11. Top: The symbols are second differences A2^, defined 
by equation (2), of low-degree (i=0,l,2) cigcnfrcquencics obtained 
from adiabatic pulsation calculations of the central model 0, and 
have the same relation to I as in Fig. 1. The solid curve is the 
diagnostic D2 determined by equations (21), (30) & (37), whose 
eleven parameters have been adjusted to fit the data by least 
squares. The measure \ 2 (mean squared differences) of the over- 
all misfit is (53nHz) 2 . The dashed curve represents the smooth 
contribution (last term in equation (21)). Bottom: Individual con- 
tributions of the oscillatory seismic diagnostic. The solid curve 
displays the He II contribution, the dotted curve is the He I con- 
tribution and the dot-dashed curve is the contribution from the 
base of the convection zone. 

within the He I ionization zone (in the case of the Sun and 
similar stars, this is not the case of He II ionization) , and the 
evanescence of the eigenfunction above the turning point 
must be taken into account. This can be achieved via the 
usual Airy-function representation. But for the purpose of 
evaluating the integral S^K, it is adequate simply to use the 
appropriate high-|i/>| sinusoidal or exponential asymptotic 
representations either side of the turning point to estimate 
the 'oscillatory' component of the integrand, which amounts 
to setting 

3 

[(dive) 2 ]osc - - cos 2t/; , for r > r t , (34) 

2 r ypcr z n 

with tp given by equation (27), without the subscripts II, 
and k by equation (26). Despite the vanishing of k at 
t — r t = (m + l)/<^, expression (34) is finite at r = r t 
because 2tp(r t ) = it/2. One can treat the evanescent region 
similarly, avoiding the singularity and making the represen- 
tation continuous at r = r t by writing 




-2 L ... 1 .... 1 .... 1 .... 1 .... 1 .... 1 .- 

1000 1500 2000 2500 3000 3500 4000 
v (fHz) 

Figure 12. Top: The symbols (with error bars computed under 
the assumption that the raw frequency errors are independent) 
represent second differences A2V, defined by equation (2), of low- 
degree solar frequencies with £=0,1,2 and 3, obtained from Bi- 
SON (Basuctal. 2006). The effective overall error in the data is 
(<r) =5.3 nHz. The solid curve is the diagnostic Diiy\ Qfc), which 
has been fitted to the data in a manner intended to provide an op- 
timal estimate of the eleven parameters a^. The values of some 
of these fitting parameters are: —S-y/'y\ Tu ~ 0.047, Tu — 819 s, 
An — 70 s, and the measure E of the overall misfit is 33 nHz. 
The direct measure x is 2.1; the minimum-x 2 fit of the function 
D2 to the data yields Xmin = 1-6- The dashed curve represents 
the smooth contribution (last term in equation (21)). Bottom: In- 
dividual contributions of the seismic diagnostic. The solid curve 
displays the He II contribution, the dotted curve is the He I con- 
tribution and the dot-dashed curve is the contribution from the 
base of the convection zone. 

where now 

i>(r) ~ \k\tu> - (m + 1) In ( !H±1 + \ K \\ + 1 . ( 36 ) 

\ TLU ) 4 

Although this expression has the wrong magnitude where 
—tp is large, that region makes very little contribution to 
the integral for 5 7 /C. We have confirmed numerically that 
the formula (35) provides a tolerable approximation. At any 
point the integrand for <5 7 /C is the product of an exponential 
and a slowly varying function F(t), say. Both F(t) and tp(r) 
can be Taylor expanded about r = fi = n + u>~ 1 ei (ei = 
en) up to the quadratic term, rendering the approximation 
asymptotically integrable in closed form. The correction to 
the result of assuming F(r) — F(fi) and tp = ipi is small, so 
we actually adopt just the leading term. 

The cyclic-eigenfrequency contribution to the entire he- 
lium glitch then becomes 

~ Au \y + \{m + 1)i/q] 
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Figure 13. Similar to Fig. 12, but with the GOLF data that 
are plotted in Fig.l. Top: The solid curve is the optimally fitted 
diagnostic D2{v;ak), with — ^t/tI-tu — 0.053, 7n ~ 810s, and 
An ~ 74 s. The measure E of the misfit is 49nHz, and \ =10.5; 



the minimum-x fit of D2 to the data yields Xn 



= 7.0. The 



dashed curve is the smooth contribution. Bottom: Individual con- 
tributions to the seismic diagnostic, with the same notation as 
Fig. 12. 



ji/J/sT 1 ( e -^ V ^ A "" 2 cos2Vi - l) 
+ Kn ' (c- 8 ^" A "" 2 cos 2i/m ~ l)] , (37) 

provided Tt < n, in which itn is denned in Section 4.1 and m 
is denned analogously; the phase ipi is given by ipi = ip(fi), 
where ip(fi) is given by equation (27), with fii replaced by 
rjfn, and with Ku replaced by ki = k(t\). If r t > n, we 
replace ki with Kij and cos2i/;i with 1 — exp(2i/'i — 7r/2), 
where now i/>i = V'('n) is given by equation (36). 

Some of the parameters obtained by fitting to the 
numerical second differences the generalization D2{v,Q-k) 
of the asymptotic formula (22) obtained by using expres- 
sion (37) for S-yV, and still using equation (30) for 8 c v, arc 
plotted in Fig. 10. A comparison of the fitting parameters 
and of x 2 between Figs 5, 6 and 10 shows that the improved 
seismic diagnostic D2 results in an even better fit to the adi- 
abatically computed eigenfrequencies than does the seismic 
diagnostic Do , which includes a description of only the He II 
ionization, and is also better than the seismic diagnostic D\ 
with only the improved treatment of the asymptotic phase 
function. 

A fit of the seismic diagnostic D2 to /S.2V of the adi- 
abatically computed eigenfrequencies of the central stellar 
model is shown in the top panel of Fig. 11, and a fit to A2^ 
of solar-frequency observations by BiSON (Basuetal. 2006) 



in the top panel of Fig. 12. It demonstrates how well the os- 
cillatory contributions to the second frequency differences of 
low-degree acoustic modes are approximated by our seismic 
diagnostic, at least for frequencies above 1300/iHz. The lower 
panels in both figures show the individual contributions from 
Hel Hell and from the discontinuity in d 2 lnp/dr 2 at the 
base of the convection zone. 

It seems likely that because the ionization diagnostic, 
and also the signature of the base of the convection zone, 
are greatest at the lowest frequencies, it is important for 
our purpose that the frequencies of the gravest modes of the 
star be measured; indeed, those modes are the least influ- 
enced by the convective fluctuations near the stellar surface, 
and therefore their frequencies are in principle the most pre- 
cisely defined, although they are apt to have low amplitudes 
and consequently be difficult to detect. We note that the 
GOLF data do extend to frequencies lower than the others, 
and with smaller standard errors. These modes dominate 
the fitting, as is evident from a comparison of the effective 
overall error (a) of the modes in the entire GOLF data set 
and that of the subset within the MDI frequency range (see 
Section 2); the higher frequencies are relatively poorly de- 
termined. Consequently in the error-weighted fits illustrated 
in Figs 1 and 13, the gross properties of both the ionization 
signature and the signature of the base of the convection 
zone are determined predominantly by the lowest-frequency 
modes. However, we must emphasize that the asymptotic 
representation used to interpret the seismic diagnostic be- 
comes less reliable at lower frequencies, and that therefore 
any model calibration might be more severely contaminated 
by structural properties other than that caused directly by 
He- induced ionization. 

It would be unwise to attempt to infer Y directly from 
the values of the fitting parameters obtained from this in- 
vestigation because the reference models were not computed 
with the most up-to-date microphysics; however, they are 
adequate for providing the partial derivatives needed for the 
calibration of a more realistic model. 

One might wonder what the effect would be of includ- 
ing the degree dependence of the mode frequencies on the 
oscillatory component of Sv. It can be taken into account by 
using the full leading-order asymptotic expression for the 
vertical component K of the wavenumber, which is given by 
equation (B2), and making the corresponding corrections to 
the eigenfunctions. For the purposes of estimating the nu- 
merators in equations (6) and (B9) it is adequate to ignore 
the buoyancy frequency and expand K just to leading order 
in I. More care must be taken in estimating the measures 
X and X* of the inertia, because the integrals extend down 
to the lower turning points of the modes where the degree- 
dependent term in K is as important as the other terms. 
The overall effect on 8 osc v is quite small, although, as we 
discuss in the next section, the i-dependence is not entirely 
negligible. Indeed, by including the degree-dependence in 
the fitting procedure it is possible to reduce x 2 further, but 
only by about 10 per cent. The improvement, as can be seen 
in a plot analogous to Fig. 11, is noticeable only at low fre- 
quencies, for which the lower turning points are higher in the 
star and have a larger impact on the modes. However, the 
resulting modifications to the parameters <57/7| TII , Tn, An 
and t c are tiny. Therefore we pursue this matter no further. 

The results of this investigation indicate that the 
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asymptotic expression for A2V can be calibrated against ob- 
servation to obtain tolerable estimates of aspects of the effect 
on 7 of helium ionization. They therefore augur the success 
of a future asteroseismic calibration. 



5 DISCUSSION 

Although we have decided to ignore the degree-dependence 
in our fitting, it is at least instructive to enquire how the 
degree I influences the ionization signature in the second 
differences. To this end we study a model in which we omit 
the effects of the base of the convection zone, accounting for 
contributions to the second frequency differences from only 
the ionization zones. Such conditions are represented by an 
adiabatically stratified stellar model satisfying the following 
equations: 

dp Gmp 
dr r 2 ' 



dm . 2 

— — = VKT p , 

dr 



dp 
dr 



Gmp 2 
7pr 2 



(38) 
(39) 
(40) 



where m is the mass enclosed in the sphere of radius r (not 
to be confused with the polytropic index which it repre- 
sents elsewhere in this paper) . In constructing the model we 
adopted 7 = 7(p,p, Y, Z) obtained from the EFF equation of 
state (Eggleton, Faulkner & Flannery 1973). The usual reg- 
ularity condition m/r 2 ~ ► as r — > and the boundary con- 
dition p = p c were adopted at the centre, and the conditions 
m = M,p — p s were imposed at the surface r = R, where 
p a is obtained from model S of Christensen-Dalsgaard et al. 
(1996). There results an eigenvalue problem for the eigen- 
value p c . 

First, we concentrate on the smooth contributions to 
the second frequency differences. We represented them by a 
third-degree polynomial in Sections 2 and 4, as in the last 
term in equation (21). The symbols in the upper left panel 
of Fig. 14 show second differences A2V of adiabatically com- 
puted eigenfrequencies for a helium-free adiabatically strat- 
ified reference model (Y re f = 0, Z = 0.02). The dotted lines 
connect symbols of like degree I. They differ only at the low- 
est frequencies, and by a relatively small amount (values for 
I — 4 lie the furthest from a single curve, as expected, but 
modes of such high degree are not considered elsewhere in 
this paper). In the upper right panel are plotted the sec- 
ond differences A25v of the oscillatory component associ- 
ated with helium ionization of a model in which the helium 
abundance Y is 0.25 (Z — 0.02). They were obtained from 
evaluating numerically the formulae (6)-(8) obtained from 
the variational principle, with £7 = 7|y=o.25 — 7|y=o.o and 
with the numerically computed eigenfunctions £ (the depth- 
dependence of ^7, plotted as a function of acoustic depth t, 
is illustrated in the lower right panel of Fig. 14), and are 
essentially independent of degree I. With the help of equa- 
tion (5) the second frequency differences of the Y = 0.25 
model can be obtained by taking the sum of A2V of the 
Y = model and A28V. The result is depicted in the lower 
left panel of Fig. 14. Once again the values associated with 
I — 0,1,2 and 3 fall almost on a single curve. Therefore 
ignoring I is quite a good procedure, as is indeed evident 



from the fit to the second differences of the frequencies of 
the solar model presented in Fig. 11. That fit is not per- 
fect, however, and it may be that it could be improved by 
extracting the /-dependence from the smooth contribution 
using the asymptotic relation (20) . The Z-dependence of the 
oscillatory component is weaker; it arises principally from 
the Z-dependence of the inertia T. 

On the whole, the amplitudes <57/7|m of the acoustic 
glitch inferred by fitting the asymptotically derived formulae 
to the frequency differences A2V of the test models described 
in Section 3.1 agree reasonably well with the actual values 
in the model, especially when our preferred expression for 
the phase function is used. The amplitude increases with Y, 
as it should, even when the He ionization zones are repre- 
sented by only a single Gaussian function (at least for the 
model sequence in which age varies). We notice in passing, 
however, that Monteiro & Thompson (2005) report that if 
the triangular approximation is made to the depression in 
7, then the inferred amplitude varies in the wrong sense. 
The other parameters in our fitting are not reproduced so 
well as the amplitude. The estimated widths An are rather 
small, which is not what one would expect at least from 
Fig. 8, because the actual perturbation to 7 appears to be 
somewhat broader than the Gaussian approximation used 
to represent it, although not by as much as the differences 
suggested by comparing Figs 9 and 10; it is comforting, how- 
ever, that when the two helium ionizations are treated sepa- 
rately, as they should be, the magnitude of the discrepancy 
in An is the least. Also in error are the acoustic depths of 
the Hell ionization zones and the bases of the convection 
zones, although they are better when the improved phase 
function is employed, as one would expect. We suspect that 
the remaining errors with the improved phase function have 
arisen because the actual models, and real stars too, arc 
not polytropes: also there is some ambiguity in obtaining 
an appropriate location of the effective acoustic surface of 
the Sun, and hence the origin of the coordinate r, because 
that depends on precisely how the layers in the convection 
zone near the upper turning points influence the dynamics 
of the modes, and consequently how the region adopted for 
the extrapolation of c 2 or w 2 should be chosen; moreover, 
the deviation of the stratification from the polytropic be- 
haviour influences the relation between tp and our polytropic 
approximation throughout the convection zone, resulting in 
the errors in ipn and ip c being different. Finally, we note 
that the value of A c inferred from the seismic calibration 
is 20% too low, which is puzzling because we would expect 
the asymptotic expressions used to relate the discontinuity 
in the second density derivative to the amplitude of the os- 
cillatory signature to work rather well. Perhaps a significant 
source of the discrepancy is our neglect of the influence of the 
discontinuity in oj a on the structure of the oscillation eigen- 
functions (see Appendix B). In any case, we stress that the 
precision of the asteroseismic calibrations that are planned 
does not rely directly on the precision of the asymptotic ex- 
pressions. Those expressions have been used solely to design 
seismic signatures of certain properties of the stratification 
of a star to be as free from contamination by other prop- 
erties as can reasonably be expected; their role is to enable 
one to interpret the results of the calibrations, and to help 
in the assessment of the accuracy of the inferences. 

It is noteworthy that the slopes of the solid and dashed 
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Figure 14. Calculated second frequency differences for adiabatically stratified models. The upper left panel shows the results for a 
reference model with helium-free gas (Y r cf = 0.0, Z = 0.02). The upper right panel shows the second differences of the frequency 
perturbation Su according to equations (6)-(8) {Su = Slu/2tt) with £7 = 7|y=0.25 — 7|y=0.0 being the difference in 7 between the 
reference model computed with Y IC { = 0.0 and a model with Y = 0.25; £7 is depicted in the lower right panel. The sum of A2 u and 
A2 Su is plotted in the lower left panel; it is an estimate of the second frequency differences of the Y = 0.25 model. The meanings of the 
symbols in the first three panels are indicated in the first panel. All results are obtained with the numerically computed cigenfunctions 
of the (unperturbed) reference model. 



lines representing the behaviour of the two different se- 
quences of solar models in the top three panels of Figs 4, 
5, 6, 9 and 10 are not parallel. The variations of m and An 
are both small, and can be influenced by details of the up- 
per superadiabatic boundary layer in the convection zone, 
for example, which is not of principal concern here. But one 
might naively expect, and indeed hope, the values of 5~//~/\t u 
to depend almost solely on the value of Y, for then the inter- 
pretation would have been simple. But that is evidently not 
the case. A substantial direct influence of heavy elements 
on the depression of 7 by Hell ionization may be partially 
responsible. But it is likely that the predominant influence 
comes about through the dependence on the specific entropy 
s of the adiabat deep in the convection zone, which is deter- 
mined by the (^-dependent) calibration of the models to the 
solar radius and luminosity. We reiterate that nevertheless 
this result does not compromise the precision of the astero- 
seismic calibration, just its interpretation. It is no surprise 
that the slopes of A C (Y) and t c (Y) are different for the two 
sequences of solar models, because they depend on the in- 
terfacing of the convection zone with the radiative interior, 
and the latter is quite sensitive to the opacity, and hence 
to the value of Z (and, therefore in a calibrated sequence of 
solar models, to Y). 

Our investigation has been carried out with a single 
equation of state. Although it can hardly be doubted that 
our signature of the magnitude of the depression in 7 is quite 



robust, it is certainly the case that relating that to Y is less 
so. The value of S~//~/\ ru is determined locally by the values 
of Y and s, and the magnitude of the error in how they are 
related depends directly on the error in the equation of state 
at tu, although it is also due partly to the inaccuracy of our 
representation of the glitch (see Fig. 7). The error would 
certainly be adequately small for the asteroseismic calibra- 
tions of the foreseeable future. But to determine the value 
of Y does seem to require knowledge of s, which can be de- 
termined only by an additional diagnostic. That diagnostic 
might have been dependent on uncertainties in conditions 
far from the Hell ionization zone, which may be difficult to 
assess. But fortunately the immediate environment of the 
helium ionization zone is isentropic, so one needs to deter- 
mine only a constant, rather than a function: given s (and 
the gravitational acceleration, which hardly varies through 
the helium ionization zones of solar-like stars, and to a high 
degree of precision may be considered to be constant too), 
the functional form of 7 is rather well constrained, and a sim- 
ply measurable property of it, such as An, might be used 
to obtain another, different, relation between Y and s. That 
would enable each of these two quantities to be determined. 
The constancy of s disengages the calibration from uncer- 
tainties elsewhere in the star. (We expect the value of tii to 
be less useful for these purposes, for it depends on conditions 
in the vicinity of the upper turning points of the modes; it 
should, however, be useful as a diagnostic of the treatment of 
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convection, which determines the structure of the superadi- 
abatic boundary layers.) However, before such a calibration 
is carried out, more thorough tests of the reliability of the in- 
ferred value of An must be performed. Essentially the same 
suggestion has been made already by Monteiro & Thompson 
(2005) to distinguish between equations of state, although 
they too offered no cogent reason to suggest whether it would 
provide a reliable test. If the hypothesis that we are putting 
forward is correct, this measure of the specific entropy would 
not alone necessarily provide a robust criterion for selecting 
any particular equation of state. However, it would certainly 
limit the uncertainty in the 5^-Y relation needed for the as- 
teroseismic calibration. 

We should also point out that the variation of 7 in the 
He ionization zone would be influenced by the presence of a 
magnetic field, which we have ignored in this study. Evidence 
for low-amplitude solar-cycle variations has been seen in the 
raw frequencies of a combination of solar models of both low 
and intermediate degree (e.g. Goldreich et al. 1991), which 
relate directly to conditions in the Hell ionization zone 
(Gough 1994), and, more recently, in fourth frequency dif- 
ferences (Basu & Mandel 2004, 2006). Such variations could 
be greater in other stars, and would therefore noticeably 
contaminate the low-degree diagnostic. Unless adequately 
accounted for, they would degrade any seismic calibration 
of Y. 

We emphasize that our aim in this work is not to in- 
vent a functional representation of the second frequency dif- 
ferences, the diagnostic of our choice, or of any other diag- 
nostic for that matter, with the intention of merely being 
able to fit the data well. Instead, it is principally to ob- 
tain a formula that can be related directly to properties of 
the stratification of the star that we believe to be perti- 
nent to the abundance of helium. Although the degree to 
which we succeed in fitting the data provides some measure 
of our success, it is not the sole criterion guiding our choice 
of function. Indeed, functions that are not obviously directly 
related to the stratification but which are tailored to mimic 
the data should be almost bound to fit the data better. But 
evidently they must be less reliable for diagnosis. It is en- 
couraging, however, that our function actually does fit the 
frequency data better than any other that has been tried 
before (Houdek & Gough 2006), particularly when the data 
are fitted over a large range of frequency, although we has- 
ten to add that its ability to measure the acoustic glitches 
is not perfect. In this regard it is interesting to note that 
the fit to the GOLF data of D2, illustrated in Fig. 13, is no 
better than the fit of D\, illustrated in Fig. 1, yet the values 
inferred for —5^h\ Tll , m and An are likely to be superior. 

Finally, we emphasize that although our analysis is 
based on asymptotic theory, the actual calibrations would 
be carried out using numerically computed eigenfrequencies, 
and their precision would be independent of the accuracy of 
the asymptotics. As is often the case in asteroseismology, 
the role of the asymptotics is to motivate the design of a 
calibration procedure, and to facilitate its interpretation in 
physical terms; it is not used in the calibration itself. 



6 SUMMARY AND CONCLUSION 

We have designed asteroseismic signatures of helium ioniza- 
tion which we trust will provide a reliable basis for calibrat- 
ing stars against a suitable grid of theoretical models, using 
second frequency differences A2V n ,l = fn-1,1 — 2f„,( + fcVj+i,i- 
Those signatures can be used to determine the magnitudes, 
widths and locations of the depressions £7 of the first adia- 
batic exponent 7 due to He ionization. 

The degree to which our original simple diagnostic (22) 
appears to reproduce the magnitude and width of the de- 
pression in 7 in the He II ionization zone can be judged by 
comparing the inferred values plotted in Fig. 5 with those 
measured from the models in Fig. 4. The values of the ampli- 
tudes of the depression are overestimated by some 50%, and 
the widths underestimated by a similar amount, yielding an 
approximately correct magnitude of the glitch integral, rn- 
However, the slopes of the dependence of those values on he- 
lium abundance Y, particularly amongst models with vary- 
ing heavy-element abundance Z, are not reproduced. In this 
regard it should be appreciated that the reason may actually 
be that our procedure for characterizing the amplitudes and 
widths of the depressions in the models may not be seismi- 
cally appropriate. The variation of the acoustic depth rn of 
the centre of ionization is more-or-less reproduced. However, 
the actual magnitude appears to be discrepant; that is due in 
part to our crude representation of acoustic phase, and per- 
haps in part to an uncertainty in the origin of r. Roughly 
half of that discrepancy is accounted for by our improved 
treatment of the phase, discussed in Section 4.1, as can be 
seen in Fig. (6). That improvement also brings the amplitude 
of the 7 depression almost into agreement with the model 
values, but the improvement in the inferred widths An is 
only slight. We suspected that the problem with the widths 
arises from our explicit neglect of He I ionization in our di- 
agnostic formula, in which we represented the two distinct 
helium depressions, which in reality separately influence the 
oscillation eigenfrequencies, by a single Gaussian function. 
However, we note that the variation of An is less than the 
variation of the other ionization-related quantities, and may 
therefore perhaps not be a robust diagnostic. Nevertheless, 
we have introduced a separate representation of He I ioniza- 
tion in Section 4.2, the result of which can be assessed by 
comparing with model values in Fig. 9 with the seismically 
inferred counterparts plotted in Fig. 10. Interestingly, the 
slope with respect to Y of the measured model values of rn 
for the Z- varying sequence are changed, whereas that of the 
inferred values is not. Moreover, the slope of the amplitudes 
of the inferred 7 depressions in the Z-varying sequence is 
improved only slightly. But the magnitude of the discrep- 
ancy between the inferred and the actual values of An is 
reduced by a factor three. 

We have found that, for solar-like stars at least, the 
magnitude of ($7 increases approximately linearly with he- 
lium abundance Y, a property which has assisted us in ac- 
counting for the He I ionization zone by relating it directly 
to the contribution from Hell ionization. By so doing, a 
somewhat more faithful overall representation of the Hell 
ionization is achieved by the seismic calibration. The mag- 
nitude of the depression in 7 is not a function of Y alone, but 
depends also on at least one more parameter. We argue that 
for a solar calibration there is just one, namely the specific 
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entropy; for other stars it must depend also on the gravi- 
tational acceleration, whose determination would be part of 
the overall calibration. If the situation is indeed that sim- 
ple, then a potential asteroseismic calibration procedure is 
almost in hand. 
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APPENDIX A: SIGNATURE OF RAPID 
VARIATION OF THE ADIABATIC EXPONENT 

The adiabatic eigenfrequencies w of a nonrotating star in the 
Cowling approximation satisfy a variational principle which 
can be written in the form (Chandrasekhar 1963) 

- 2 = % (Al) 



where 



and 



47T 



pt-tdV, 



(A3) 



the integrals being over the volume of the star. Omitted 
from equation (Al) is a possible surface integral; that inte- 
gral contributes to lo at most a slowly varying function (with 
respect to lo), which in any case in practice is small, so we 
have ignored it at the outset for simplicity. By adopting the 



Cowling approximation we have also ignored the contribu- 
tion from the gravitational-potential perturbation, because 
that too is only slowly varying. For low-degree modes of high 
order, the integrand in K, is dominated by the first term 
(whose integral we call K,\), except possibly in localized re- 
gions where the equilibrium density p varies abruptly, such 
as at the base of the convection zone. Therefore the contri- 
bution to Slo from a region of relatively rapid variation of 7 
produced by the ionization of an abundant element may be 
approximated by 



Slo ■ 



SK-i — lo 2 5T 



5Ki 
2lo1 



2lo1 

since \6Ki\ > lo 2 \51\, where 



(A4) 
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= : 5yJCi + SplCi 



and 
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SpS-tidV , 



(A6) 



(A7) 



in which £7 and Sp are to be regarded as the differences be- 
tween 7 and p in the equilibrium star and those in a similar 
smooth model in which ionization is absent. In view of the 
variational principle, the differences (presumed to be small) 
between the values of £ in the star and in the smooth model 
do not contribute to Slo to leading order in small quantities. 
If 57 is defined to be nonvanishing only in the ionization 
region where it varies rapidly with r, then one expects <5 7 /Ci 
to contribute the dominant oscillatory (with respect to fre- 
quency) component to lo, because Sp, which is related to £7 
via its radial derivative, varies with r more slowly, as does 
Sp which in the helium ionization zone is directly related 
to Sp through the adiabatic constraint. (Note that for low- 
amplitude perturbations Sp , Sp, 8-y, the relation between Sp 
and Sp does not depend on £7 to leading order.) We have 
confirmed numerically that this is indeed the case. Conse- 
quently, most of the oscillatory contribution to Slo, in the 
absence of rotation, is contained in the term {2loI)~ 1 8 1 K,\. 

In estimating the integrals X and K, it is adequate for 
our purposes to use the high-order, low-degree asymptotic 
expressions for the eigenfunctions of the smooth model, and 
then to evaluate the integrals asymptotically. In solar-like 
stars the He II ionization zone is well inside the propagating 
region of the acoustic modes, and the leading term in the 
Liouville- Green approximation to the solution of the adia- 
batic wave equation suffices, as is described in Section 2. 
That is not necessarily the case for the He I ionization zone, 
and formally one should retain the Airy-function form of the 
JWKB approximation. However, for the purpose of evaluat- 
ing K,, we have found it sufficient to adopt the large-phase 
sinusoidal or exponential representations either side of the 
upper turning point r t of the mode, even up to the turning 
point itself, but scaling the eigenfunction div£ in the evanes- 
cent zone with a constant chosen to render that eigenfunc- 
tion continuous at r — r t . 

Rotation influences the eigenfrequencies via the Corio- 
lis force (with respect to a rigid frame of reference rotating 
locally with the fluid) and via centrifugal effects which both 
influence the oscillation dynamics directly and distort the 
equilibrium structure of the star from sphericity. One can 
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estimate the influence with an analysis similar to that pre- 
sented here using the more general variational principle of 
Lynden-Bell & Ostriker (1967). In particular, one would ob- 
tain an estimate of the oscillatory contribution to modes of 
like degree I and azimuthal order m induced by a tachocline 
flow near the base of the convection zone of an essentially 
axisymmetric star. However, if lo is regarded as the multi- 
plet frequency, obtained as the equally weighted average of 
the frequencies of all the modes of like I and n, but varying 
azimuthal order, such contributions are suppressed, at least 
when rotation is slow (e.g. Gough 1993). Then the multiplet 
frequencies provide a good diagnostic of helium ionization. 
However, if the centrifugal force is not negligible, its average 
effect on the equilibrium structure of the star would need to 
be taken into account in the grid of reference stellar models 
used for comparison. 



APPENDIX B: SIGNATURE OF THE 
ACOUSTIC GLITCH AT THE BASE OF THE 
CONVECTION ZONE 

The acoustic glitch at the base of the convection zone is es- 
sentially a discontinuity in the second derivative of the tem- 
perature, associated with which are discontinuities in the 
second derivatives of sound speed and density. That renders 
all three terms in the integrand for K,, given by equation 
(A2), susceptible to perturbations of the same order. Then 
terms must be calculated with some care, because there can 
be cancellation in leading order when the terms are com- 
bined (e.g. Monteiro, Christensen-Dalsgaard & Thompson 
1994). It is more straightforward, however, to work di- 
rectly (in the planar Cowling approximation) with the single 
second-order differential equation for the r-dependent fac- 
tor p(r) in the Lagrangian pressure perturbation eigenfunc- 
tion, for then only one quantity is discontinuous, namely 
the acoustic cutoff frequency co a (and, of course, the sec- 
ond derivative of the eigenfunction) . That equation is (e.g. 
Gough 1993) 



d 2 ^ 9 



(Bl) 



where ~ rp p. The vertical wavenumber is given by 
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K 2 = LO -LO*. 
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in which N is the buoyancy (Brunt- Vaisala) frequency: 
1 9 



N"=g 



H 



also 



2 

= 



4H 2 



J dr 



c 

4JP 



c 2 d 2 In p 
2 dr 2 



(B3) 



(B4) 



where g is the local acceleration due to gravity and H = 
— (dlnp/dr) -1 is the density scale height. 

Equation (Bl) is to be solved subject to a regularity 
condition at r = and an appropriate causality condition 
at r = R. It is written in self-adjoint form, and therefore 
it follows immediately that for modes with frequencies well 
below the acoustic cutoff frequency characteristic of the at- 



mosphere, the equation 
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(B5) 



defines a variational principle for lo 2 amongst acceptable 
functions ^; those functions are such that the boundary 
terms that arise from the integration by parts to obtain the 
first term in equation (B5) either vanish (at r = 0) or can 
be neglected (at r — R) because they contribute at most a 
slowly varying function of lo. At the base of the convection 
zone, r — r c , lo' 2 suffers a discontinuity 



r 2i ' 



(B6) 
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Let lo 2 be extrapolated smoothly from the convection 
zone down into the radiative interior to define lo 2 , with the 



property w a 

r 2 2 



lo' as r ■ 



0. Additionally, let us define 



(B7) 



which satisfies 5u), = for r > r c , Slo„ 



-A c at r 



and Suig, — > as r — > 0. Granted that <5w a is small enough for 
linearization to be valid, the contribution S c lo 2 to Slo 2 from 
the discontinuity in lo 2 is then given by 

1) f R 8N 2 
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Formally it was necessary to include the second and third 
terms on the left-hand side of this equation because the vari- 
ational principle holds only for perturbations in the model 
that preserve hydrostatic support, and changing lo 2 leads 
to a change 8N 2 in iV 2 . However, these terms are small for 
low-degree high-order p modes (and zero for radial modes). 
Consequently they can safely be ignored, yielding 



S c LO 



Io'SloI 



c- 2 * 2 dr 



/ c- 2 * 2 dr 



8JC 



(B9) 



This formula is adequate for the purposes of studying the 
signature of helium ionization. We note, however, that the 
linearization may not be adequate for studying the base of 
the convection zone; for that endeavour one may need to 
account for the perturbations to the eigenfunctions, as did 
Roxburgh and Vorontsov (1994). From asymptotic theory 
the high-order eigenfunctions can be approximated by (e.g. 
Gough 1993): 

*o^ _1/2 sin^, (BIO) 



between the turning points rb and r t at which K, vanishes, 
where 



tp = / if dr + - ; 



(Bll) 



then 
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Figure Bl. Acoustic cutoff frequency w a in the vicinity of the 
base of the convection zone. The solid curves show ui 2 of the 
central model and the dashed curve is the sum of the solid curve 
and of equation (B15) evaluated with tq = 80s. 
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in which r t = r(r t ) and r b = r(r h ) are the acoustic depths 
of the upper and lower turning points. In our introductory 
analysis in Section 2 we ignored cj a ; in that case the last of 
equations B12 follows immediately from the equation above. 
When o) a is included in the polytropic approximation, as in 
Section 4.1, the first-order correction to vanishes, and 
therefore equation (B12) still holds. We have also ignored 
the i-dependence of K, in the penultimate approximation in 
equations (B12), consistently with our reduction of equation 
(B8) to equation (B9). At this level of approximation, 
becomes the same as the leading-order approximation (13) 
to the inertia X defined by equation (8) after setting >]/o = w- 



x cos [2-0C + tan 1 (2tooj)] 



(B17) 



Finally, we notice that only the second term in equation 
(B4) is discontinuous, whence 



A c = 
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"~2 



d 2 lnp 



dr 2 



(B18) 



at the base of the convection zone, 



Setting (cK) 2 i 
finally obtain for the oscillatory contribution 5 c ,osc<^ to ui 
from the discontinuity in the acoustic cutoff frequency at 
r = r c : 





-1/2 
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x cos [2ip c + tan -1 (2Too>)] , 



(B19) 



where c c = c(r c ). 

Within our initial approximation discussed in Section 2, 
in which the acoustic cutoff frequency cj a was not explicitly 
taken into account, we set ip c = t c u) + e c , where e c is con- 
stant, as in equation (17). When w a is taken into account 
as in Section 4, the phase is tp c = tp{f c ), in which tp(f c ) 
is essentially given by equation (27) with fii replaced by 



4>c 



T c +u> e c : 

, -,\ _ i / m + 1 \ 7T 
t K C T C LJ — (m + 1) cos — ; I + — • 
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(B20) 



where k c = [1 — {m + l) 2 /f^u 2 ] 1 ^ 2 , which is tantamount 
to adding a frequency-dependent phase term to the con- 
stant e c . For stars with acoustically deep convection zones, 
the correction to the initial approximation is small, and the 
frequency dependence is therefore weak. In the case of the 
sun, for example, f c cu ~ 40 for a 3 mHz mode, whence 
dipc/duj — K C T C ~ 0.995f c . 



The integral 8JC in equation (B9) can be approximated by 

8JC ~ *o / Suj^c" A" _1 sin i/)dr 

~ i*o / " Su;l(cKy 2 [1 - cos 2^] K dr. (B13) 

Once again it is adequate to set K ~ lo/c, and to set V = 
i> c + ui(t — t c ) with 



f rt TT 



(B14) 



Furthermore, beneath the convection zone we approximate 
lo 2 by an exponential function: 



-A c e 



-(t-t c )/t 



(B15) 



where To <S T is determined by fitting expression (B15) to 
the central model of our grid; it is compared with the model 
in Fig. Bl. Then, since 5u> 2 varies much more rapidly than 
(cK)~ 2 - in fact cK ~ u) = constant - one may ignore 
the variation of the latter over a few to. It follows that the 
oscillatory component of S C IC is then 

f e -0— tc)/tc cos2 ^dr (B16) 

1/2 
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APPENDIX C: SECOND DIFFERENCE OF THE 
OSCILLATORY FREQUENCY COMPONENT 



Consider an oscillatory signal of the form 

Su n = A n COSXn 



(CI) 



with v n = L0 n /2TT, where A n — A n (uj n ) and where x n — 
x n (u! n ) varies almost linearly with u) n throughout the range 
of frequencies u>„ that we are considering. For clarity we have 
suppressed the index I. To estimate the second frequency 
difference A2<5^ n = 5v n +i — 28v n + 8v„-\ we regard ui n as 
a continuous function of n and expand x n ±i and A n ±i in 
Taylor series about x n and A n respectively, retaining only 
two terms for x n ±i and three for A n ±i: 

x n±1 ~x n ±a, A n±1 ~ (l±a + b)A n , n > 1 , (C2) 

where uio = 2-kuq approximates the large (angular) fre- 
quency separation and 
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d 2 A ri 



The derivatives of u) n are obtained by differentiating equa- 
tion (20) with respect to n at fixed I; it is consistent with 
the rest of this analysis to retain only the leading term: 
dLO n /dn ~ loo, d 2 uj„/dn 2 ~ 0. Substituting equations (C2) 
into the definition of the second difference yields 

A2<5i/„ = [(f + a + b) cos(x n + a) — 2 cos:r n 
+ (1 — a + b) cos(x n — a)]A n 
= FA n cos(x„ — 3) , 



where 

F = 2{[1 - (1 + b) cos a} 2 + a 2 sin 2 a} 1 ' 2 
and 

5 = tan" 



f — (f + b) cos Q 



(C5) 
(C6) 

(C7) 



Had the frequency separation Ai^ n been small compared 
with the scale of variation on A n , then A\8un and A2<5^n 
would approximate the first and second derivatives of 5v n 
with respect to n; but it is not, and the formulae (C5)-(C7) 
differ substantially from the second derivative. 

In the first approximation to the contribution from 
Hell ionization given by equation (15), A„ oc Lj n e~ 2Au " n 
and x n = 2(rno; n + en). Therefore, 



a ~ 2u>otii 
and 



(C8) 



a ~ -(4A 2 I o; 2 1 - Tjuouj^ . 
b ~ 2A 2 I w 2 (4A 2 I u; 2 -3). 



(C9) 

In the case of the Sun, lo ~ 8.8 x 10 _4 s _1 , An ~ 80s and 
Tii ~ 800 s; whence a ~ 1.4. In the middle of the frequency 
range used for our calibration of solar models, tu n ~ 27r x 2.7 
mHz ~ 1.7 x l(T 2 s _1 ; whence a ~ -0.33 and b ~ 0.04. The 
latter is so small partly because of a near cancellation of the 
two terms in parentheses at this frequency (cancellation is 
exact at v ~ 2 mHz), but it never exceeds this value by more 
than a factor four in the frequency range of interest, and is 
greatest at the highest frequency where the amplitude A n is 
the smallest. The effect of the neglected higher-order terms 
in the Taylor expansions of A n and x n is even smaller. The 
value of the amplitude factor F ~ 1.7, and S ~ —0.4. 

When the acoustic cutoff and the second-order term 
in the asymptotic eigenfrequency relation (20) are taken 
into account, as in Section 4.1 , the amplitude factor for 
the ionization glitch becomes A„ oc AuK^[u> n + (m + 
l)wo/2] exp(-2« n A n u;£) and x„ = 2ip u , where Vn = ip(m), 
with fii = Tn + uj~ €u and ip(fn) given by equation (27). 
Then 



a ~ 2wotiikii ; 

Ku = /t(fn) is given by equation (26), and 



-[4A n w 2 + (t n 2 
2A? I wg[4A? I w^- 



- 1 - (1 + /3 a ) ]W0W„ 
1-2(1 + /3a)- 1 ], 



(C10) 



(Cll) 



where /3 a = \{m + ljwow^ 1 . We have simplified the ex- 
pression for & by ignoring a term (m + l) 2 /4A n f I 2 I c 1 j* , and 
other yet smaller terms, which typically augment the al- 
ready small b by only about 1 per cent. For the He I con- 
tribution, one must replace Ku and An by ki and Ai when 



Table CI. Relative signal error magnification due to order-fc 
differencing of frequencies, evaluated at v=l.b, 2.75 and 4.0 mHz. 

ek/Fk 



He II ionization zone 
=1.5 2.75 4.0 mHz 



Base of convection zone 
u=1.5 2.75 4.0 mHz 



1 


1.38 


1.30 


1.33 


0.83 


0.81 


0.80 


2 


1.91 


1.46 


1.30 


0.71 


0.73 


0.74 


3 


3.22 


2.44 


2.33 


0.76 


0.77 


0.76 


4 


4.74 


2.90 


2.27 


0.70 


0.75 


0.76 


5 


7.98 


4.92 


3.90 


0.79 


0.81 


0.81 


6 


12.10 


5.79 


3.65 


0.74 


0.81 


0.83 


8 


31.43 


11.19 


5.57 


0.81 


0.90 


0.93 


10 


51.36 


21.95 


9.27 


0.89 


1.02 


1.05 



Tt < n. When n is in the evanescent zone of the mode, 
so that 6vn = A n exp(— x n ), then A25t> n = FA n exp(— x n ), 
where F = 2[(1 + b) cosh a — 1 — a sinh a] . The derivatives a 
and 6 are still given by equations (C10) and (Cll), but now 
a ~ 2woft|Ki|- 

At the base of the convection zone, x n = 2ip c + 
tan - (2row n ) — 2tp c and a ~ 2ojot c because k c = k(t c ) 



is very close to unity. Here, f c — r c + ui 
A n oc u;- 2 (l + l/4r 2 cj 2 )- 1/2 , whence 



2o; ^ 1 (l + l/8r 2 o; 2 )/(l + l/4r 2 o. 2 ) 
2 (3 + 5/8r 2 w 2 )/(l + l/4r 2 o. 2 ) 2 . 
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Moreover, 



(C12) 



Again we have simplified b, this time by ignoring 1/48tooj* 
in the parentheses in the numerator, which typically con- 
tributes about 0.2 per cent. In the case of the Sun, to ~ 80 s 
and r c ~ 2273s (see Fig. Bl). Whence a ~ -0.10 and 
b ~ 0.01 in the middle of the frequency range, and F ~ 3.4 
and S ~ 0.05. 

It is interesting to estimate the signal-to-noise ratio for 
differences AkSvn of various orders k. It is straightforward 
to extend the analysis above to obtain formula of the kind: 
AkSfn = FkA n cos(xn — 5k)- The frequency differences in- 
volve a weighted sum of fc + 1 frequency increments 8v„, 
which may be written 

fe/2 

A k 5v„ = ^2 dklSVn+l , (C13) 
!=-fc/2 

where, for example, d 2 i = (1, —2, 1) and d 4; = 
(1,-4,6,-4,1). Presuming, for simplicity, that the statis- 
tical uncertainties in the frequency measurements are inde- 
pendent, all with standard deviation a, then the variance 
of the possible error in Ak8v n is e|a 2 = d|;(T 2 . Hence 
the magnification of the measurement error, in units of the 
magnification of the magnitude of the oscillatory signal, is 
ek/F k . That quotient depends on the nature and location 
of the acoustic glitch and the central frequency of the con- 
tributing modes, through Fk, and, through dki, on the order 
fc of the difference. It is tabulated in Table CI for the Hell 
ionization zone and for the base of the convection zone, at 
the median frequency v n = 2.75 mHz and the two extreme 
frequencies 1.5 and 4.0 mHz. For the Hell ionization zone 
the signal amplitude factor Fk is augmented by a factor of 
about 1.4 when k is increased by unity, but that is more than 
offset by the augmentation, by a factor of about 1.9, of the 
error magnification . Choosing which order fc of the differ- 
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ence to adopt must therefore require a tradeoff between the 
value of the error quotient et /Fk and the degree of difficulty 
in separating the smooth and oscillatory components. It is 
interesting to note that at the base of the convection zone 
€k/Fk is lowest for k = 2 (except, marginally, at the lowest 
frequency), although not greatly so, which perhaps favours 
A2<5^ as a diagnostic, in accord with the conclusion of Ballot 
et al. (2004) drawn from Monte Carlo calculations. 



APPENDIX D: EVALUATION OF THE 
CONVECTION-ZONE DISCONTINUITY 

The difference equations that were used to represent the 
differential equations of stellar evolution (which are based 
simply on second-order-accuracy centred differences in space 
and time on a fixed Lagrangian spatial mesh) do not ad- 
equately account for the mathematical properties of the 
structure in the vicinity of the base of the convection zone. 
The convective heat flux was computed from a local mixing- 
length formalism with a mixing length set to a constant pro- 
portion of a pressure scale height (which therefore does not 
vanish at the edge of the convection zone), which leads to 
a discontinuity in the second spatial derivative of the den- 
sity, and hence, according to equation (B4) , in the critical 
cutoff frequency, which, as is common, was implicitly pre- 
sumed to be absent when adopting the difference equations 
for the computations. Therefore particular care should be 
taken when estimating conditions at the base of the con- 
vection zone for the purpose of determining the parame- 
ters r c and A c in the seismic signature. Perhaps the most 
straightforward way to have computed the structure would 
have been to ensure that there was always a mesh point at 
the base of the convection zone and to tailor the difference 
scheme either side of that point to accommodate the correct 
mathematical structure of the solutions. But the models had 
already been computed, and so it was necessary instead to 
fit appropriate functional forms to the solutions close, but 
not too close, to the base of the convection zone for the 
purpose of evaluating the magnitude of the discontinuity. 

The functional forms adopted were those derived by 
Christensen-Dalsgaard, Gough & Thompson (1991). Above 
the base of the convection zone, r = r c , the density varia- 
tion is smooth, and a polynomial extrapolation of dlnp/dr 
and d 2 lnp/dr 2 , the latter of which was calculated by nu- 
merical differentiation, to r — r c is adequate. For r < r c we 
fitted to dlnp/dlnr the function 
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(2m - h)y 
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in which z = p/po and 
R 



h p c GM 
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withp c =p(r c ), pc = p(r c ); also / = (4 + A + v)/(3 + A), h = 
f-l,q = Xhw/{1 + A/), y = (1 + 2A)(/ + \w)hw/[2(2Xf + 
f + I)], u = (1 + A)/, w — 1/(3 + v), po being an integration 
constant to be determined in the fitting process, where for- 
mally A = (d In x/d hip),} and v = — (d In x/d lni9) p are the 
logarithmic partial derivatives of opacity x with respect to 
density p and temperature both derivatives being evalu- 
ated at r — r c . Expression (Dl) was fitted by least squares 




Figure Dl. Density derivatives near the base of the convection 
zone. The top panel shows the first derivative of density. The plus 
symbols display the density derivatives of the original (unmodi- 
fied) model and are connected by the dotted curve. The vertical 
dotted line indicates the location of the base of the convection 
zone of the unmodified model. The dashed curve is the result of 
fitting the analytical expression Dl by least squares to the plus 
symbols of the original (unmodified) model within the fitting do- 
main indicated by the horizontal solid line. The square symbols 
show the density derivatives of the modified model and are con- 
nected by the solid curve (they lie on the dashed curve above the 
mid-point of the fitting domain, and on the dotted curve below 
it). The vertical solid line indicates the new location of the base 
of the convection zone, r c . The bottom panel shows the result for 
the second density derivatives of the modified model obtained by 
numerical differentiation. 

to the solar models in the region indicated in Fig. Dl by ad- 
justing the parameters A, po, A, and v, the outcome of which 
is illustrated as the dashed curve in the upper panel of the 
figure for the central model 0. Its intersection with the ex- 
trapolation of din p/dlnr in the convection zone determines 
r c , from which it can then be confirmed that equation (D2) 
is satisfied. From this can be computed the second deriva- 
tive of In p, which is illustrated in the lower panel of Fig. Dl, 
and hence the discontinuity d 2 In p/dr 2 | rc + — d 2 In p/dr 2 | rc _ 
which is needed for evaluating A c according to equation (18). 
Both A c and r c = r(r c ) of the nine test models are plotted 
against Y in the first lower two panels of Fig. 4. 



